; This file reads the line intensity tables and the gadget output of a
; simulation and calculates the emission from particles. It then creates
; emission maps and saves them in HDF5 format.
;
;------------------------------------------------------------------------------

   pro map, sim=sim, snap=snap, line=line, nslices=nslices, $
            save=save, pixel=pixel, thick=thick, zconst=zconst, $
	    cut_temp=cut_temp, mintemp=mintemp, maxtemp=maxtemp, $
	    cut_den=cut_den, minden=minden, maxden=maxden, $
	    zoom=zoom, npix=npix, delta=delta, xc=xc, yc=yc, zin=zin

;------------------------------------------------------------------------------
;
;   By Serena Bertone   2007
;
;------------------------------------------------------------------------------

   print, ' usage: map'
   print
   print, ' Keywords:'
   print, '    sim :     simulation to use	    (REF_L100N128, default)'
   print, '    snap:     desired snapshot number    (snap=30, default)'
   print, '    line:     desired line	            (line=O 6 1032A, default)'
   print, '    pixel:    pixel size in arcsec	    (15 arcsec, default)'
   print, '    thick:    thickness of slice	    (20 Mpc/h, default)'
   print, '    save :    saves postscript output    (PS not saved, default)'
   print, '    nslices:  number of slices	    (1, default)'
   print, '    zconst :  use constant metallicities (no, default)'
   print, '    zin:      initial z coord. for slice (0.0, default)'
   print, '    zoom:     zoom in on selected region (all slice, default)'
   print, '    npix:     number of pixels for zoom  (250, default)'
   print, '    delta:    width of region for zoom   (0.1, default)'
   print, '    xc:       x central coord. for zoom  (0.3606, default)'
   print, '    yc:       y central coord. for zoom  (0.42, default)'
   print, '    cut_temp: consider only particles with given temperatures'
   print, '    mintemp:  minimum temperature in Log (5, default)'
   print, '    maxtemp:  maximum temperature in Log (7, default)'
   print, '    cut_den:  consider only particles with given density'
   print, '    minden:   minimum density in Log	    (3, default)'
   print, '    maxden:   maximum density in Log	    (4, default)'
   print
   print, ' mintemp and maxtemp work only when keyword cut_temp is set!'
   print, ' minden and maxden work only when keyword cut_den is set!'
   print

   if (not keyword_set(sim))      then sim     = 'REF_L100N128'
   if (not keyword_set(snap))     then snap    = 30
   if (not keyword_set(line))     then line    = 'O  6  1032A'
   if (not keyword_set(pixel))    then pixel   = 15.0
   if (not keyword_set(thick))    then thick   = 20.0
   if (not keyword_set(nslices))  then nslices = 1
   if (not keyword_set(save))     then save    = 0 else save = 1
   if (not keyword_set(zconst))   then zconst  = 0 else zconst = 1
   if (not keyword_set(zin))      then zin     = 0.0
   if (not keyword_set(zoom))     then zoom    = 0 else zoom = 1
   if (not keyword_set(npix))     then npix    = 250
   if (not keyword_set(delta))    then delta   = 0.1
   if (not keyword_set(xc))       then xc = 0.3606
   if (not keyword_set(yc))       then yc = 0.42
   if (not keyword_set(cut_temp)) then cut_temp = 0 else cut_temp = 1
   if (not keyword_set(cut_den))  then cut_den = 0 else cut_den = 1
   if (cut_temp eq 1 and cut_den eq 1) then cut_temp = 0
   if (cut_temp eq 1 or cut_den eq 1)  then cut = 1 else cut = 0
   if (not keyword_set(mintemp) and cut_temp eq 1) then mintemp = 5.0
   if (not keyword_set(maxtemp) and cut_temp eq 1) then maxtemp = 7.0
   if (not keyword_set(minden) and cut_den eq 1) then minden = 3.0
   if (not keyword_set(maxden) and cut_den eq 1) then maxden = 4.0
   if (save eq 0)   then print, ' Map not saved as postscript.' 
   if (zconst eq 1) then print, ' Using constant metallicities.' 
   if (zoom eq 0)   then print, ' Zooming in selected region: xc=',xc,' yc=',yc
   print, ' Simulation: ', sim
   print, ' Snapshot:   ', snap
   print, ' Line:       ', line
   if (cut_den eq 1) then begin
      print, '    minden: ', minden
      print, '    maxden: ', maxden
   endif
   if (cut_temp eq 1) then begin
      print, '    mintemp: ', mintemp
      print, '    maxtemp: ', maxtemp
   endif
   print

   ; Directories with input files
   ; Gadget files:
   mem = memory()
   dir = '/disks/strw17/owls/sims/'
   snap = strcompress(snap, /remove_all)
   file = dir + sim + '/data/snapshot_0' + snap + '/snap_0' + snap + '.0.hdf5'
   ; Emissivity tables
   dir_tables = '/disks/strw17/serena/Tables/'
   ; List of redshifts of the tables
   redshift_file = '/home/bertone/Cloudy/Data/UVtables/redshifts.dat'

   ; Read the gadget files
   print, ' Read the snapshots...'
   ; Read some attributes
   file_id = H5F_OPEN(file)
      group_id = h5g_open(file_id,"/Header")
	 hubble = read_attribute_hdf5(group_id,"HubbleParam")
	 redshift = read_attribute_hdf5(group_id,"Redshift")
	 omega = read_attribute_hdf5(group_id,"Omega0")
	 omegabaryon = read_attribute_hdf5(group_id,"OmegaBaryon")
	 omegalambda = read_attribute_hdf5(group_id,"OmegaLambda")
	 boxsize = read_attribute_hdf5(group_id,"BoxSize")
	 numpart_total = read_attribute_hdf5(group_id, "NumPart_Total")
	 numfiles = read_attribute_hdf5(group_id, "NumFilesPerSnapshot")
      H5G_close,group_id
      group_id = h5g_open(file_id,"/Constants")
	 protonmass = read_attribute_hdf5(group_id,"PROTONMASS")
	 solarmass  = read_attribute_hdf5(group_id,"SOLAR_MASS")
	 PI         = read_attribute_hdf5(group_id,"PI")
	 cm_per_mpc = read_attribute_hdf5(group_id,"CM_PER_MPC")
	 gravity    = read_attribute_hdf5(group_id,"GRAVITY")
	 lightspeed = read_attribute_hdf5(group_id,"C")
	 planck     = read_attribute_hdf5(group_id,"PLANCK")
      H5G_close,group_id
      group_id = h5g_open(file_id,"/Parameters/ChemicalElements")
	 CalciumOverSilicon = read_attribute_hdf5(group_id,"CalciumOverSilicon")
	 SulphurOverSilicon = read_attribute_hdf5(group_id,"SulphurOverSilicon")
      H5G_close,group_id
   H5F_CLOSE, file_id
   ngas = numpart_total(0)
   ; Check that nslices and thick are consistent
   if (nslices gt 1) then $
      if (nslices * thick gt boxsize) then nslices = fix(boxsize / thick)
   print, ' thickness of slice:	', thick, ' Mpc/h'
   print, ' number of slices:	', nslices
   print, ' hubble:   ', hubble
   print, ' redshift: ', redshift
   print, ' number of particles: ', ngas
   print

   ; Calculate luminosity distance
   if (redshift lt 1.e-2) then rcom = boxsize / (2.0 * hubble) $
   else rcom = comovingdistance(redshift, omega, omegalambda, hubble)
   lum_dist = rcom * (1.0 + redshift) * cm_per_mpc
   lum_dist = alog10(4.0 * PI) + 2.0 * alog10(lum_dist)

   ; Calculate dimension of box in arcsec
   angdiam = rcom / (1.0 + redshift)              ; Mpc
   pixel_arcsec = temporary(pixel)
   pixel_rad = pixel_arcsec * PI / (180.0 * 3600.0)
   theta_box_rad = boxsize / (hubble * angdiam * (1.0 + redshift))
   theta_box_arcsec = theta_box_rad * 180.0 * 3600.0 / PI
   npixel = fix(theta_box_arcsec / pixel_arcsec)
   omega_box = 2.0 * PI * (1.0 - cos(theta_box_rad / 2.0))
   omega_pixel = 2.0 * PI * (1.0 - cos(pixel_rad / 2.0))
   print, ' pixel dimension: ', pixel_arcsec, ' arcsec'
   print, ' pixel dimension: ', pixel_rad, ' rad'
   print, ' pixel dimension: ', omega_pixel, ' sr'
   print, ' box dimension: ', theta_box_arcsec, ' arcsec'
   print, ' box dimension: ', temporary(theta_box_arcsec)/3600.0, ' deg'
   print, ' box dimension: ', theta_box_rad, ' rad'
   print, ' box dimension: ', omega_box, ' sr'
   print, ' npixel: ', npixel
   if (zoom eq 1) then begin
      theta_box_rad = boxsize * delta / (hubble * angdiam * (1.0 + redshift))
      pixel_rad = theta_box_rad / npix
      pixel_arcsec = pixel_rad * 180.0 * 3600.0 / PI
      omega_pixel = 2.0 * PI * (1.0 - cos(pixel_rad / 2.0))
      npixel = npix
      print
      print, ' pixel dimension: ', pixel_arcsec, ' arcsec'
      print, ' pixel dimension: ', pixel_rad, ' rad'
      print, ' pixel dimension: ', omega_pixel, ' sr'
      print, ' selected box dimension: ', theta_box_rad, ' rad'
      print, ' npixel: ', npixel
   endif
   print

   ; Read the particle data...
   print, ' Read the particle data...'
   temperature = read_gadget_hdf5(file, 0, 'Temperature', silent=1)
   density     = read_gadget_hdf5(file, 0, 'Density', silent=1)
   mass        = read_gadget_hdf5(file, 0, 'Mass', silent=1)
   hydrogen    = read_gadget_hdf5(file, 0, 'Hydrogen', smooth=1, silent=1)
   sfr         = read_gadget_hdf5(file, 0, 'StarFormationRate', silent=1)

   ; Read element abundance
   el = strmid(line, 0, 2)
   if (el eq 'He') then element = 'Helium'
   if (el eq 'C ') then element = 'Carbon'
   if (el eq 'N ') then element = 'Nitrogen'
   if (el eq 'O ') then element = 'Oxygen'
   if (el eq 'Ne') then element = 'Neon'
   if (el eq 'Mg') then element = 'Magnesium'
   if (el eq 'Fe') then element = 'Iron'
   if (el eq 'Si' or el eq 'Ca' or el eq 'S ') then element = 'Silicon'
   if (zconst eq 0) then begin
      abund = read_gadget_hdf5(file, 0, element, smooth=1, silent=1)
      w = where(abund eq 0)
      if (w(0) ne -1) then abund(w) = 1.0e-15
      w = 0
      ; Calcium and Sulphur are calculated relative to the silicon abundance.
      if (el eq 'Ca') then abund = temporary(abund) * SulphurOverSilicon
      if (el eq 'S ') then abund = temporary(abund) * CalciumOverSilicon
      print, ' done.'
      print
   endif

   ; Read the emissivity table
   ; Emissivity tables are expressed as arrays with indexes:
   ; (line_lambda, hden, temp)
   print, ' Read emissivity tables...'
   z_tables = fltarr(49)
   openr, unit, redshift_file, /get_lun
      readf, unit, z_tables
   free_lun,unit
   w = value_locate(z_tables,redshift)
   if (redshift-z_tables(w) le z_tables(w+1)-redshift) then $
      ztable = temporary(z_tables(w)) $
   else ztable = temporary(z_tables(w+1))
   ztable = strcompress(ztable, /remove_all)
   dir_tables = dir_tables + strmid(ztable, 0, 6) + '/Tables/'
   if (el eq 'He') then name_element = 'helium'
   if (el eq 'C ') then name_element = 'carbon'
   if (el eq 'N ') then name_element = 'nitrogen'
   if (el eq 'O ') then name_element = 'oxygen'
   if (el eq 'Ne') then name_element = 'neon'
   if (el eq 'Mg') then name_element = 'magnesium'
   if (el eq 'Si') then name_element = 'silicon'
   if (el eq 'Fe') then name_element = 'iron'
   if (el eq 'Ca') then name_element = 'calcium'
   if (el eq 'S ') then name_element = 'sulphur'
   table_file = dir_tables + name_element + '.hdf5'
   print, ' Read emissivity table: ', table_file
   file_id = H5F_OPEN(table_file)
      hden = read_dataset_hdf5(file_id,'logd')
      temp = read_dataset_hdf5(file_id,'logt')
      ; Element solar abundances in cloudy units.
      ; Number densities relative to hydrogen.
      ; Compilation from literature.
      group_id = H5G_open(file_id,"/header/abundances")
	 h_solar  = read_attribute_hdf5(group_id,'hydrogen_solar')
	 he_solar = read_attribute_hdf5(group_id,'helium_solar')
	 c_solar  = read_attribute_hdf5(group_id,'carbon_solar')
	 n_solar  = read_attribute_hdf5(group_id,'nitrogen_solar')
	 o_solar  = read_attribute_hdf5(group_id,'oxygen_solar')
	 ne_solar = read_attribute_hdf5(group_id,'neon_solar')
	 mg_solar = read_attribute_hdf5(group_id,'magnesium_solar')
	 si_solar = read_attribute_hdf5(group_id,'silicon_solar')
	 fe_solar = read_attribute_hdf5(group_id,'iron_solar')
	 ca_solar = read_attribute_hdf5(group_id,'calcium_solar')
	 s_solar  = read_attribute_hdf5(group_id,'sulphur_solar')
      H5G_close, group_id
      table_all = read_dataset_hdf5(file_id,'lines')
      lines = read_dataset_hdf5(file_id,'lambda')
   H5F_CLOSE, file_id

   nden  = n_elements(hden)
   ntemp = n_elements(temp)
   hdenmin = hden(0)
   hdenmax = temporary(hden(nden-1))
   tempmin = temp(0)
   tempmax = temporary(temp(ntemp-1))
   table = fltarr(nden,ntemp)
   w = where(temporary(lines) eq line)
   table(*,*) = temporary(table_all(w,*,*))
   w = 0
   table = reform(temporary(table))
   print, ' Table read.'

   ; Atomic weights (atomic mass) - taken from the periodic table 
   m_h  = 1.00794
   m_he = 4.002602
   m_c  = 12.0107
   m_n  = 14.0067
   m_o  = 15.9994
   m_ne = 20.1797
   m_mg = 24.3050
   m_si = 28.0855
   m_fe = 55.845
   m_ca = 40.078
   m_s  = 32.065

   ; Make cuts on particle numbers
   tt = alog10(temporary(temperature))
   if (cut eq 1) then begin
      if (cut_temp eq 1) then $
         wcut = where(tt ge mintemp and tt le maxtemp)
      if (cut_den eq 1) then begin
	 meanrhob = 3.0 * omegabaryon * hubble^2 / (cm_per_mpc * 1e-7)^2 / $
                   (8.0 * PI * gravity) * (1.0 + redshift)^3
	 meanrhob = float(meanrhob)                       ; proper g/cm^3
	 den = alog10(density / meanrhob)
	 wcut = where(den ge minden and temporary(den) le maxden)
      endif
      ngas     = n_elements(wcut)
      tt       = temporary(tt(wcut))
      sfr      = temporary(sfr(wcut))
      mass     = temporary(mass(wcut))
      abund    = temporary(abund(wcut))
      density  = temporary(density(wcut))
      hydrogen = temporary(hydrogen(wcut))
      print, ' number of particles left after cuts: ', ngas
   endif

   ; Prepare for the interpolation
   emis  = fltarr(ngas)
   lhden = alog10(density * hydrogen / protonmass)
   ; Factor for conversion
   factor = 2.0 * lhden + alog10(temporary(mass) * m_h / (density $
            * temporary(hydrogen)))
   ; Remove particles with SFR>0 and density higher than a threshold
   wsfr = where(temporary(sfr) gt 0.0 or lhden gt -1.0)

   if (el eq 'He') then cc = m_he * he_solar
   if (el eq 'C ') then cc = m_c  * c_solar
   if (el eq 'N ') then cc = m_n  * n_solar
   if (el eq 'O ') then cc = m_o  * o_solar
   if (el eq 'Ne') then cc = m_ne * ne_solar
   if (el eq 'Mg') then cc = m_mg * mg_solar
   if (el eq 'Si') then cc = m_si * si_solar
   if (el eq 'Fe') then cc = m_fe * fe_solar
   if (el eq 'Ca') then cc = m_ca * ca_solar
   if (el eq 'S ') then cc = m_s  * s_solar

   ; Interpolate
   print, ' interpolate the emissivity tables...'
   emis = interpolate(temporary(table), (temporary(lhden) - hdenmin) * $
	  (nden - 1.0) / (hdenmax - hdenmin), (temporary(tt) - tempmin) * $
	  (ntemp - 1.0) / (tempmax - tempmin))
   emis = reform(temporary(emis))
   if (zconst eq 0) then $
   emis = temporary(emis) + temporary(factor) + alog10(temporary(abund) / cc) $
   else emis = temporary(emis) + temporary(factor) - 1.0

   ; Transform emission in erg/(s*cm^2)
   emis = temporary(emis) - lum_dist

   ; Convert to photons by dividing by the energy per photon.
   lwave = strlen(line)
   wave = float(strmid(line, 5, lwave-1))
   emis = temporary(emis) - alog10(planck * lightspeed / (wave * 1e-8) $
          / (1.0 + redshift))
	  ; photons/cm^2/s

   ; Star forming particles have zero emission
   if (wsfr(0) ne -1) then emis(temporary(wsfr)) = -100.0
   print, ' done.'
   unit = ' photons/cm^2/s/sr'
   print, ' total emission from box: ', total(10^emis, /nan)/omega_box, unit
   print

;nn = n_elements(emis)
;close, 1
;openw, 1, '/home/bertone/Owls/Cosmop/CPtest/myemission.dat'
;   printf, 1, nn
;   printf, 1, emis
;close, 1

;------------------------------------------------------------------------------

   ; Define plot properties
   !p.charsize = 1.2
   !p.charthick = 3
   !p.thick = 3
   !x.thick = 3
   !y.thick = 3
   !p.multi = 0
   loadct, 5

   ; Define plot region such that there is room for a color bar.
   xwindow = [0.05,0.95]
   ywindow = [0.15,0.95]
   barlength = 0.8 * (xwindow[1] - xwindow[0])
   barcenter = 0.5 * (xwindow[1] - xwindow[0]) + xwindow[0]
   barposition = [barcenter - 0.5 * barlength, 0.08, $
                 barcenter+0.5 * barlength, 0.12]
   barcharsize = 1
   ; Establish plot region.
   position = [xwindow[0], ywindow[0], xwindow[1], ywindow[1]]
   px = xwindow * !d.x_vsize         ; Size of window in device pixels
   py = ywindow * !d.y_vsize

   ; Write title of plot
   if (pixel_arcsec ge 60.0) then begin
      unit = ' arcmin'
      pixel = fix(pixel_arcsec / 60.0)
   endif else begin
      unit = ' arcsec'
      pixel = fix(pixel_arcsec)
   endelse
   pixel = strcompress(pixel, /remove_all)
   zz = redshift
   zz = strcompress(zz, /remove_all)
   zz = strmid(zz, 0, 4)
   title = line + '  z=' + zz + '  pixel=' + pixel + unit

   ; Read particle positions
   print, ' read particle positions and sph smoothing kernel...'
   allpos = read_gadget_hdf5(file, 0, 'Coordinates', silent=1)
   allkernel = read_gadget_hdf5(file, 0, 'SmoothingLength', silent=1)
   if (cut eq 1) then begin
      allpos = temporary(allpos(*,wcut))
      allkernel = temporary(allkernel(wcut))
   endif
   allpos = temporary(allpos) * hubble * (1.0 + redshift) / cm_per_mpc / boxsize
   allkernel = temporary(allkernel) * hubble * (1.0 + redshift) / $
               (cm_per_mpc * boxsize)
   print, ' done.'

;------------------------------------------------------------------------------

   ; Find directory containing cosmoplotter
   dir = '/home/bertone/Codes/Joop/cosmoplotter'
   Object_File = dir + '/HsmlAndProject_perbc.so'
   print, ' SPH executable is in: ', dir
   print, ' SPH executable is: ', Object_File

   ; For those particles that have zero smoothing length
   DesNgb = 32L

   ; Loop over all slices
   box = 1.0
   for islice=0, nslices-1 do begin
      print
      print, ' Doing slice ', islice+1, '...'
      ; Limits on z-axis
      llimit = zin / boxsize
      ulimit = (zin + thick) / boxsize

      ; Select coordinate ranges for zoom
      if (zoom eq 1) then begin
	 ; Fix the x coordinate
	 dx1 = xc - delta / 2.0
	 dx2 = xc + delta / 2.0
	 if (dx1 lt 0.0) then begin
            dx1 = 0.0
	    dx2 = delta
	 endif
	 if (dx2 gt 1.0) then begin
	    dx2 = 1.0
	    dx1 = 1.0 - delta
	 endif
	 ; Fix the y coordinate
	 dy1 = yc - delta / 2.0
	 dy2 = yc + delta / 2.0
	 if (dy1 lt 0.0) then begin
	    dy1 = 0.0
	    dy2 = delta
	 endif
	 if (dy2 gt 1.0) then begin
	    dy2 = 1.0
	    dy1 = 1.0 - delta
	 endif
      endif else begin
         dx1 = 0.0
	 dx2 = 1.0
	 dy1 = 0.0
	 dy2 = 1.0
      endelse

      ; Number of particles and pixels in the calculation:
      n0 = npixel
      w = where(allpos(2,*) ge llimit and allpos(2,*) le ulimit)
      np = n_elements(w)
      ; Generate coordinate array for particles
      xyz = fltarr(3,np)
      xyz[0,*] = allpos(0,w)
      xyz[1,*] = allpos(1,w)
      xyz[2,*] = allpos(2,w)
      xyz = temporary(xyz)
      ;Xmin = 0.0
      ;Xmax = 1.0
      ;Ymin = 0.0
      ;Ymax = 1.0
      Xmin = dx1
      Xmax = dx2
      Ymin = dy1
      Ymax = dy2
      Zmin = llimit
      Zmax = ulimit
      ; Projection along z-axis
      Axis1 = 0l
      Axis2 = 1l
      Axis3 = 2l
      ; Smoothing kernel and maximum used smoothing length
      Hsml = reform(allkernel(w))
      Hmax = max(Hsml)
      ; Emission vector
      value = reform(10^emis(w))
      weight = value
      w = 0.0

      ; Projected weights (normally projected mass map)
      ResultW = fltarr(n0,n0)
      ; Projected `Weight' weighted quantity
      ResultQ = fltarr(n0,n0)

;print, ' number of particles np: ', np
;print, axis1, axis2, axis3
;print, Xmin, Xmax, Ymin, Ymax, Zmin, Zmax
;print, Hmax
;print, max(xyz), max(Hsml), max(weight), max(value)

      print, ' Calling external routine ' + Object_File + $
             ' written by Volker Springel'
      S = CALL_EXTERNAL(Object_File, $
                	'findHsmlAndProject', $
                	long(np),       $
                	float(xyz),     $
                	float(Hsml),    $
                	float(weight),  $
                	float(value),   $
                	float(Xmin),    $
                	float(Xmax),    $
                	float(Ymin),    $
                	float(Ymax),    $
                	float(Zmin),    $
                	float(Zmax),    $
                	long(n0),       $    ; # pixels in x-direction
                	long(n0),       $    ; # pixels in y-direction
                	long(DesNgb),   $
                	long(Axis1),    $
                	long(Axis2),    $
                	long(Axis3),    $
                	float(Hmax),    $
                	1d0,            $    ; Simulation box size
                	float(ResultW), $
                	float(ResultQ))
      slice = transpose(ResultW)
      ; Free memory
      xyz = 0.
      Hsml = 0.
      weight = 0.
      value = 0.
      print, ' SPH interpolation done.'
      print,  ' total slice after interpol: ', total(slice)
      print

      ; Divide by solid angle of pixel
      print, ' total emission in selected box: ', total(slice, /nan), $
             ' photon/s/cm^2'
      slice = temporary(slice) / omega_pixel
      print, ' total emission in selected box: ', total(slice, /nan), $
             ' photon/s/cm^2/sr'
      max = max(slice)
      min_slice = 1.e-15
      min = min_slice
      zero = where(slice lt min, nzero)
      if (nzero gt 0) then slice(zero) = min
      slice = alog10(temporary(slice))
      slice = float(temporary(slice))
      print, ' minimum set to', min_slice

      ; Prepare for output
      dir = '/disks/strw17/serena/Owls/' + sim
      oscommand = 'mkdir ' + dir
      spawn, oscommand, osoutput
      iline = strcompress(line, /remove_all)
      pixel = fix(pixel_arcsec)
      pixel = strcompress(pixel, /remove_all)
      ii = islice + 1
      ii = strcompress(ii, /remove_all)
      idelta = strcompress(delta, /remove_all)
      file = 'map_'
      base = file
      if (cut_temp eq 1) then $
         base = '/tcut_' + strcompress(fix(mintemp), /remove_all) + '_' + $
	        strcompress(fix(maxtemp), /remove_all) + '_' + file
      if (cut_den eq 1) then $
         base = '/dcut_' + strcompress(fix(minden), /remove_all) + '_' + $
	        strcompress(fix(maxden), /remove_all) + '_' + file
      if (zoom eq 1) then base = 'zoom_' + idelta + '_' + file
      if (zconst eq 1) then base = 'zconst_' + file
      filehdf = dir + '/' + base + snap + '_' + iline + '_pix' + pixel + $
                '_slice' + ii + '.hdf5'

      if (save eq 1) then begin
	 ; Plot image
	 image = slice
	 ma = 2.0
	 mi = -6.0
	 image = bytscl(temporary(image), min=mi, max=ma, top=255)

	 set_plot, 'PS'
	 fileout = dir + '/' + base + snap + '_' + iline + '_pix' + pixel + $
                   '_slice' + ii + '_all.ps'
	 device, filename = fileout, xsize=20, ysize=22.5, /color, $
        	 xoffset=1, yoffset=5, bits=8, /portrait

	 ; Null ticks
	 ti = [' ', ' ', ' ']
	 xmin = 0.0
	 xmax = boxsize
	 ymin = 0.0
	 ymax = boxsize
	 plot, findgen(5), /nodata,  $
               title = title, $
	       xtickname=ti, ytickname=ti, xticks=2, yticks=2, $
	       xrange = [xmin, xmax], xstyle=1, $
	       yrange = [ymin, ymax], ystyle=1, $
	       position=position
	 tv, temporary(image), px(0), py(0), xsize=px(1)-px(0), $
	     ysize=py(1)-py(0)

	 ; Draw the scale dimension
	 nt = 8.0
	 yy = (ymax - ymin) / 10.0
	 cc = yy * 0.3
	 plots, [xmax*2.0/nt,6.0*xmax/nt], [yy,yy], linestyle=0, $
	        thick=6, color=255
	 xyouts, 3.5*xmax/nt, yy+cc, color=255, charsize=2, $
        	 strmid(strcompress(boxsize * delta / 2.0, /remove_all), 0, 2) $
		 + ' Mpc/h'

 	 ; Colour bar
	 str = ' Log!D10!N S!DB!N (photons cm!E-2!N s!E-1!N sr!E-1!N)'
	 cbar, vmin=mi, vmax=ma, cmin=0, cmax=255, color=color, $
	       position=barposition, xcharsize=barcharsize, xtitle=str
	 device, /close_file
      endif

      ; Write the slice info in a file for future use
      units = 'Emission in Log photons/s/cm2/sr, pixel in arcsec, z in Mpc/h'
      struct_pix  = {_NAME:'pixel_arcsec',_TYPE:'attribute',_DATA:pixel_arcsec}
      struct_u    = {_NAME:'units',       _TYPE:'attribute',_DATA:units}
      struct_line = {_NAME:'line',        _TYPE:'attribute',_DATA:line}
      struct_box  = {_NAME:'boxsize',     _TYPE:'attribute',_DATA:boxsize}
      struct_red  = {_NAME:'redshift',    _TYPE:'attribute',_DATA:redshift}
      struct_zmin = {_NAME:'zmin',        _TYPE:'attribute',_DATA:zmin}
      struct_zmax = {_NAME:'zmax',        _TYPE:'attribute',_DATA:zmax}
      struct_thi  = {_NAME:'thick',       _TYPE:'attribute',_DATA:thick}
      struct_min  = {_NAME:'min_slice',   _TYPE:'attribute',_DATA:min_slice}
      struct_header = {_NAME:'header',_TYPE:'group', S1:struct_u,    $
                                                     S2:struct_line, $
                                                     S3:struct_pix,  $
                                                     S4:struct_red,  $
                                                     S5:struct_zmin, $
                                                     S6:struct_zmax, $
                                                     S7:struct_thi,  $
                                                     S8:struct_box,  $
						     S9:struct_min}
      struct_image = {_NAME:'image',_TYPE:'dataset',_DATA:slice}
      H5_CREATE, filehdf, {_TYPE:'group', S1:struct_header, $
                                          S2:struct_image}

      ; Reset the minimum coordinate of the slice
      zin = zin + thick
   endfor
   ; Free memory
   slice = 0
   emis = 0
   allpos = 0
   allkernel = 0

;------------------------------------------------------------------------------

   mem = memory()
   max_memory = mem(3) / 1.0e6
   print
   print, ' ******  Memory usage statistics:  ******'
   print, ' Max memory used by code: ', max_memory, ' Mb'
   print, ' Memory in use now:       ', mem(0) / 1.0e6, ' Mb'
   print
   print, ' Finished. Bye!!!'

end
